## Get interactive session ##
# srun --time=08:00:00 --mem=40G -p int --pty bash
# module purge;source /camp/stp/babs/working/software/modulepath_new_software_tree_2018-08-13;module load pandoc/2.2.3.2-foss-2016b;ml R/3.6.0-foss-2016b-BABS;R;
# sbatch --time=08:00:00 --wrap "module purge;source /camp/stp/babs/working/software/modulepath_new_software_tree_2018-08-13;module load pandoc/2.2.3.2-foss-2016b;ml R/3.6.0-foss-2016b-BABS;Rscript runCM.r" --job-name="rCM" --mem=42G -o rCM.slurm >> commands.txt
# sbatch --time=18:00:00 --wrap "module purge;source /camp/stp/babs/working/software/modulepath_new_software_tree_2018-08-13;module load pandoc/2.2.3.2-foss-2016b;ml R/3.6.0-foss-2016b-BABS;Rscript runCM.r" --job-name="rCM" -p hmem --mem=300G -o rCN.slurm >> commands.txt
module purge;source /camp/stp/babs/working/software/modulepath_new_software_tree_2018-08-13;module load pandoc/2.2.3.2-foss-2016b;ml R/3.6.0-foss-2016b-BABS;R;
############################################################################### Perform integrated analysis ##
if (length(Obio@sampleDetailList) > 1) {
DefaultAssay(OsC_CF) <- "integrated"
} else {
Obio@parameterList$singleCellClusterString <- gsub("integrated", "RNA", Obio@parameterList$singleCellClusterString)
}
# Run the standard workflow for visualization and clustering This will scale on the most variable features only
OsC_CF <- ScaleData(OsC_CF, verbose = FALSE)
redVec <- names(OsC_CF@reductions)
if (!("pca" %in% redVec)) {
OsC_CF <- RunPCA(OsC_CF, npcs = Obio@parameterList$singleCellSeuratNpcs4PCA, verbose = FALSE)
}
# t-SNE and Clustering
## Add PCA clusters to data collection ##
if (!("umap" %in% redVec)) {
OsC_CF <- RunUMAP(OsC_CF, reduction = "pca", dims = 1:20)
}
if (!("tsne" %in% redVec)) {
OsC_CF <- RunTSNE(OsC_CF, reduction = "pca", dims = 1:20)
}
OsC_CF <- FindNeighbors(OsC_CF, reduction = "pca", dims = 1:20)
cat(paste(knit(text = chnkVec, quiet = T), collapse = "\n"))
Figure 1A: umap Clusterplot with parameter 0.2 here.
Figure 1B: umap Clusterplot with parameter 0.2 here.
Figure 1A: umap Clusterplot dendrogram with parameter 0.2 here.
Figure 2A: umap Clusterplot with parameter 0.3 here.
Figure 2B: umap Clusterplot with parameter 0.2 here.
Figure 2A: umap Clusterplot dendrogram with parameter 0.3 here.
Figure 3A: umap Clusterplot with parameter 0.5 here.
Figure 3B: umap Clusterplot with parameter 0.2 here.
Figure 3A: umap Clusterplot dendrogram with parameter 0.5 here.
Figure 4A: umap Clusterplot with parameter 0.7 here.
Figure 4B: umap Clusterplot with parameter 0.2 here.
Figure 4A: umap Clusterplot dendrogram with parameter 0.7 here.
Figure 5A: umap Clusterplot with parameter 0.9 here.
Figure 5B: umap Clusterplot with parameter 0.2 here.
Figure 5A: umap Clusterplot dendrogram with parameter 0.9 here.
Figure 6A: umap Clusterplot with parameter 1 here.
Figure 6B: umap Clusterplot with parameter 0.2 here.
Figure 6A: umap Clusterplot dendrogram with parameter 1 here.
Figure 7A: umap Clusterplot with parameter 1 here.
Figure 7B: umap Clusterplot with parameter 0.2 here.
Figure 7A: umap Clusterplot dendrogram with parameter 1 here.
Figure 8A: umap Clusterplot with parameter 1.3 here.
Figure 8B: umap Clusterplot with parameter 0.2 here.
Figure 8A: umap Clusterplot dendrogram with parameter 1.3 here.
Figure 9A: umap Clusterplot with parameter 1.5 here.
Figure 9B: umap Clusterplot with parameter 0.2 here.
Figure 9A: umap Clusterplot dendrogram with parameter 1.5 here.
Figure 10A: umap Clusterplot with parameter 1.7 here.
Figure 10B: umap Clusterplot with parameter 0.2 here.
Figure 10A: umap Clusterplot dendrogram with parameter 1.7 here.
Figure 11A: umap Clusterplot with parameter 1.9 here.
Figure 11B: umap Clusterplot with parameter 0.2 here.
Figure 11A: umap Clusterplot dendrogram with parameter 1.9 here.
Figure 12A: umap Clusterplot with parameter 2.1 here.
Figure 12B: umap Clusterplot with parameter 0.2 here.
Figure 12A: umap Clusterplot dendrogram with parameter 2.1 here.
Figure 13A: umap Clusterplot with parameter 2.3 here.
Figure 13B: umap Clusterplot with parameter 0.2 here.
Figure 13A: umap Clusterplot dendrogram with parameter 2.3 here.
Figure 14A: umap Clusterplot with parameter 2.5 here.
Figure 14B: umap Clusterplot with parameter 0.2 here.
Figure 14A: umap Clusterplot dendrogram with parameter 2.5 here.
library(AUCell)
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")
dfHeatmapGenes <- read.delim(
Obio@parameterList$catRefFile,
header = T,
sep = "\t",
stringsAsFactors = F
)
geneSets <- list()
for (i in 1:ncol(dfHeatmapGenes)){
genes <- as.vector(dfHeatmapGenes[2:nrow(dfHeatmapGenes),i])
genes <- genes[genes %in% rownames(x = OsC_CF@assays$RNA)]
geneSets[[names(dfHeatmapGenes)[i]]] <- genes
}
Obio@parameterList[["cat2DplotList"]] <- geneSets
###############################################################################
## Get backdrop
exprMatrix <- as.matrix(OsC_CF@assays$RNA@counts)
#logMat <- log10(exprMatrix+1)
# When using a Seurat object #
logMat <- data.frame(OsC_CF[["RNA"]]@data)
## Load tSNE coordinates ##
cellsTsne <- data.frame(OsC_CF@reductions$umap@cell.embeddings)
## done
FNbase <- paste0("CatScatter_Rankings", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
cells_rankings <- AUCell_buildRankings(exprMatrix)
min 1% 5% 10% 50% 100%
1037.00 1150.45 1648.00 1899.50 3361.00 8212.00
dev.off()
png 2
geneSets <- Obio@parameterList$cat2DplotList
cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings, aucMaxRank=nrow(cells_rankings)*0.05)
## Select thresholds ##
FNbase <- paste0("CatScatterHist", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
set.seed(123)
cells_assignment <- AUCell_exploreThresholds(
cells_AUC,
plotHist=TRUE,
nCores=1,
assign=TRUE
)
dev.off()
png 2
## Add data to dfExpr ##
## Plot CatScatters ##
for (i in 1:length(Obio@parameterList$cat2DplotList)){
HMname <- names(Obio@parameterList$cat2DplotList)[i]
tag <- gsub("[.]", "_", HMname)
FNbase <- paste0("CatScatterHist_", HMname, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
selectedThresholds <- cells_assignment[[i]]$aucThr$thresholds
if ("minimumDens" %in% rownames(selectedThresholds)) {
pThr <- selectedThresholds["minimumDens", "threshold"]
} else if ("Global_k1" %in% rownames(selectedThresholds)){
pThr <- selectedThresholds["Global_k1", "threshold"]
} else {
pThr <- selectedThresholds[1, "threshold"]
}
if (nrow(cellsTsne) > 15000){
cex = 0.25
} else if (nrow(cellsTsne) > 1000){
cex = 0.5
} else {
cex = 1
}
## Get AUC matrix ##
tSNE.df <- data.frame(cellsTsne, cell=rownames(cellsTsne))
mAUC <- getAUC(cells_AUC)[HMname,rownames(tSNE.df)]
dfAUC <- data.frame(mAUC)
dfAUC[["cellID"]] <- row.names(dfAUC)
dfAUC <- merge(dfAUC, tSNE.df, by.x = "cellID", by.y = "cell")
input <- list(
"x_axis" = "UMAP1",
"y_axis" = "UMAP2",
"gene" = HMname
)
dotsize <- cex
legendNote <- paste0(
" The following genes of this dataset are represented in this figure: ",
paste0(sort(Obio@parameterList$cat2DplotList[[i]]), collapse = ", ")
)
plotList[[tag]] <- ggplot(data = dfAUC, aes(x=UMAP_1, y=UMAP_2, color = mAUC)
)+ geom_point( shape=16, size = dotsize
) + scale_color_gradient(low="grey", high="darkblue"
) + xlab(input$x_axis) + ylab(input$y_axis) + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12)
)+ ggtitle(paste0("Category: ", input$gene)
) + coord_fixed(ratio = 1) #+ theme(legend.position="none")
FNbase <- paste0("CatScatter", HMname, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
## Create R markdown chunk ##
figLegend <- paste0(
"**Figure ",
figureCount,
":** Category Scatter showing gene category ",
HMname, ". ", legendNote,
". Download a pdf of this figure [here](", FNrel,"). "
)
figureCount <- figureCount + 1
NewChnk <- paste0(
"### Category Feature Plot ",HMname,
"\n```{r CatFeatPlot1_",
i,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
}
save(OsC_CF,
file = paste0(
Obio@parameterList$localWorkDir,
Obio@parameterList$project_id,
".SeuratCF"
)
)
if (length(plotList) > 3){
tabVar <- ".tabset .tabset-fade .tabset-dropdown"
} else {
tabVar <- ".tabset .tabset-fade .tabset-pills"
}
cat(paste(knit(text = chnkVec, quiet = T), collapse = "\n"))
Figure 15: Category Scatter showing gene category Enteric_Glia. The following genes of this dataset are represented in this figure: S100b, Sox10. Download a pdf of this figure here.
Figure 16: Category Scatter showing gene category Enteric_Neuron. The following genes of this dataset are represented in this figure: Chat, Elavl3, Elavl4, Gad1, Nos1, Ret. Download a pdf of this figure here.
Figure 17: Category Scatter showing gene category Proliferation. The following genes of this dataset are represented in this figure: Mki67, Top2a. Download a pdf of this figure here.
Figure 18: Category Scatter showing gene category Apoptosis. The following genes of this dataset are represented in this figure: Acin1, Acvr1c, Aifm3, Akt1, Apc, Bbc3, Bcap31, Birc2, Blcap, Bmx, Bnip1, Capn10, Casp2, Casp3, Casp6, Casp7, Casp8, Cdk5rap3, Cecr2, Cflar, Cideb, Cidec, Clspn, Dedd2, Dffa, Dffb, Dicer1, Dnase1l3, Dnase2a, Dnm1l, Endog, Ern2, Fnta, H1f0, Hmgb2, Htra2, Kpna1, Kpnb1, Madd, Prkcd, Prkcq, Ptk2, Rock1, Satb1, Sharpin, Stk24, Taok1, Tardbp, Top2a, Xkr8. Download a pdf of this figure here.
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS: /camp/apps/misc/stp/babs/manual/software/r/R-3.6.0-foss-2016b/lib64/R/lib/libRblas.so
## LAPACK: /camp/apps/misc/stp/babs/manual/software/r/R-3.6.0-foss-2016b/lib64/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_GB.utf-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.utf-8 LC_COLLATE=en_GB.utf-8
## [5] LC_MONETARY=en_GB.utf-8 LC_MESSAGES=en_GB.utf-8
## [7] LC_PAPER=en_GB.utf-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.utf-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] AUCell_1.6.1 ggtree_1.16.6
## [3] sctransform_0.2.1 KernSmooth_2.23-15
## [5] fields_10.3 maps_3.3.0
## [7] spam_2.5-1 dotCall64_1.0-0
## [9] DoubletFinder_2.0.3 scales_1.1.1
## [11] RColorBrewer_1.1-2 mixtools_1.2.0
## [13] RMySQL_0.10.20 DBI_1.1.0
## [15] DT_0.14 Seurat_3.1.5
## [17] knitr_1.29 forcats_0.5.0
## [19] stringr_1.4.0 purrr_0.3.4
## [21] readr_1.3.1 tidyr_1.1.0
## [23] tibble_3.0.1 tidyverse_1.3.0
## [25] ggplot2_3.3.2 dplyr_1.0.0
## [27] DESeq2_1.24.0 SummarizedExperiment_1.14.1
## [29] DelayedArray_0.10.0 BiocParallel_1.18.1
## [31] matrixStats_0.56.0 Biobase_2.44.0
## [33] GenomicRanges_1.36.1 GenomeInfoDb_1.20.0
## [35] IRanges_2.18.3 S4Vectors_0.22.1
## [37] BiocGenerics_0.30.0
##
## loaded via a namespace (and not attached):
## [1] R.utils_2.9.2 reticulate_1.16 tidyselect_1.1.0
## [4] RSQLite_2.2.0 AnnotationDbi_1.46.1 htmlwidgets_1.5.1
## [7] Rtsne_0.15 munsell_0.5.0 codetools_0.2-16
## [10] ica_1.0-2 future_1.17.0 withr_2.2.0
## [13] colorspace_1.4-1 highr_0.8 rstudioapi_0.11
## [16] ROCR_1.0-11 listenv_0.8.0 labeling_0.3
## [19] GenomeInfoDbData_1.2.1 bit64_0.9-7 farver_2.0.3
## [22] vctrs_0.3.1 treeio_1.8.2 generics_0.0.2
## [25] xfun_0.15 R6_2.4.1 rsvd_1.0.3
## [28] locfit_1.5-9.4 bitops_1.0-6 assertthat_0.2.1
## [31] promises_1.1.1 nnet_7.3-12 gtable_0.3.0
## [34] globals_0.12.5 rlang_0.4.6 genefilter_1.66.0
## [37] splines_3.6.0 lazyeval_0.2.2 acepack_1.4.1
## [40] broom_0.5.6 checkmate_2.0.0 BiocManager_1.30.10
## [43] yaml_2.2.1 reshape2_1.4.4 modelr_0.1.8
## [46] backports_1.1.8 httpuv_1.5.4 Hmisc_4.4-0
## [49] tools_3.6.0 ellipsis_0.3.1 ggridges_0.5.2
## [52] Rcpp_1.0.4.6 plyr_1.8.6 base64enc_0.1-3
## [55] zlibbioc_1.30.0 RCurl_1.98-1.2 rpart_4.1-15
## [58] pbapply_1.4-2 cowplot_1.0.0 zoo_1.8-8
## [61] haven_2.3.1 ggrepel_0.8.2 cluster_2.0.8
## [64] fs_1.4.2 magrittr_1.5 data.table_1.12.8
## [67] RSpectra_0.16-0 lmtest_0.9-37 reprex_0.3.0
## [70] RANN_2.6.1 fitdistrplus_1.1-1 hms_0.5.3
## [73] patchwork_1.0.1 mime_0.9 evaluate_0.14
## [76] xtable_1.8-4 XML_3.99-0.3 jpeg_0.1-8.1
## [79] readxl_1.3.1 gridExtra_2.3 compiler_3.6.0
## [82] crayon_1.3.4 R.oo_1.23.0 htmltools_0.5.0
## [85] segmented_1.2-0 later_1.1.0.1 Formula_1.2-3
## [88] geneplotter_1.62.0 lubridate_1.7.9 formatR_1.7
## [91] dbplyr_1.4.4 MASS_7.3-51.4 rappdirs_0.3.1
## [94] Matrix_1.2-18 cli_2.0.2 R.methodsS3_1.8.0
## [97] igraph_1.2.5 pkgconfig_2.0.3 rvcheck_0.1.8
## [100] foreign_0.8-71 plotly_4.9.2.1 xml2_1.3.2
## [103] annotate_1.62.0 XVector_0.24.0 rvest_0.3.5
## [106] digest_0.6.25 RcppAnnoy_0.0.16 tsne_0.1-3
## [109] graph_1.62.0 rmarkdown_2.3 cellranger_1.1.0
## [112] leiden_0.3.3 tidytree_0.3.3 htmlTable_2.0.0
## [115] uwot_0.1.8 GSEABase_1.46.0 kernlab_0.9-29
## [118] shiny_1.5.0 lifecycle_0.2.0 nlme_3.1-139
## [121] jsonlite_1.7.0 viridisLite_0.3.0 fansi_0.4.1
## [124] pillar_1.4.4 lattice_0.20-38 fastmap_1.0.1
## [127] httr_1.4.1 survival_3.2-3 glue_1.4.1
## [130] png_0.1-7 bit_1.1-15.2 stringi_1.4.6
## [133] blob_1.2.1 latticeExtra_0.6-29 memoise_1.1.0
## [136] irlba_2.3.3 future.apply_1.6.0 ape_5.4